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A MULTI SPECTRAL DATA SIMULATION TECHNIQUE* 
Marwan J. Muasher and Philip H. Swain 


For remote sensing data analysis , several assumptions are 
commonly made. These assumptions are usually that the data are 
class-conditionally distributed multivariate normal and that the 
data used to train the classifier are representative of the area 
of interest. This second assumption actually has several parts. 
The assumption is made that in the process of training, all 
classes present in the scene are found, and all spectral sub- 
classes of each class are also represented in the training data. 
Furthermore, the parameters of the distribution of each subclass 
are also assumed to be known from the training data. Each pixel 
is assumed to come from one of the training classes, and also is 
assumed to be entirely of one cover type. 

In actual practice, these assumptions are not met. The number 
of spectral classes in the area is not known and clustering or 
some other method is used to determine the number of subclasses, 
in addition to estimating the statistics of those subclasses. 

Some of these methods also lead to non-normal subclasses. In 
particular, the clustering algorithm available through LARSYS 
truncates the tails of the subclass distributions and so leads 
to non-normal distributions. 

There are also questions relating to a single picture element. 
A single pixel in Landsat data covers an area approximately 80 
meters by 50 meters. More than one cover type may be present in 
this area and result in a "mixture pixel" observation. It is not 
clear how the distribution of the spectral response of mixture 
pixels can be related to the distribution of the spectral response 
of "pure pixels." 
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There has been much speculation in the remote sensing commu- 
nity as to the effect of the non-satisfaction of the basic assump- 
tions. Whenever new algorithms are brought forth, the old ques- 
tions are raised again, indicating that there is insufficient 
understanding of the interaction of the real attributes of the 
data and the theory of the algorithms. At times it is not clear 
whether a particular result is due to aspects of the algorithm 
or to the extent the data set deviates from the assumptions. 

In testing new algorithms, deviations from the assumptions 
may obscure the action of the new process. One way to clarify 
the situation is to apply the algorithm first to a data set 
satisfying the assumptions. 

Such a data set could be obtained artificially, through 
simulation. The analyst could then know: how many classes 

exist in the data; the true distributions of the classes, in- 
cluding normality if desired; the observations could really be 
independent; and no pixel would be a "mixture pixel." New algo- 
rithms could be studied on such a data set with the knowledge 
that any "strange" effects are indeed algorithm rather than data 
problems . 

In many cases where simulated data have been used in the 
past, the data were too artificial, in the sense that all 
aspects of the image were controlled, removing the natural 
variation in object size, position, and relationship which 
occur in real data. This limited the use of the simulated 
data sets in testing new algorithms. 

The natural spatial information occurring in multispectral 
data could be retained in a simulated image by spatially basing 
the simulation on a classification. It would be even better to 
base the simulated data on a digitized "ground truth" map if the 
spectral characteristics of the cover types were known. By 
basing the simulation on a classification, the number of classes, 
their exact distributions, and the class of each pixel in the 


area are known. If the classification was sufficiently accurate, 
then the spatial information held in the classification map will 
be close to the actual cover type map and actual spatial content 
of the original data. For each pixel in the area, a random 
vector distributed according to the pixel's class statistics 
could be generated. This becomes the simulated data vector. 

Statistical Background 

From the classification chosen as a basis for the simulation, 
the following are known: the number cf classes K, the set of 

classes {wi, i=l, . . . ,K) , the class distributions {f (<ui) ,i=l, . . . ,K) , 
their means and covariances {y^ and E^, i=l, . . . ,K} , the number of 
channels p, and the class of every pixel in the scene. 

From classifical statistics: 

(1) Let X:pxl, A:pxp and b:pxl. 

If X 'v N {0,I p ) , then Y = AX + b ^ N (b, AI p A T = AA T ) 
(where I p is the identity matrix having dimensionality p) . 

(2) Let E be a symmetric, positive definite matrix. Then there 

exists A, such that 

AA T - Z (A is denoted l **) 

To simulate a pixel which was a member of class i in the base 
classification, N(0,I p ) (the random vector for each pixel is 
indepedent of other vectors) is generated. (See Appendix I.) 

Next Y = E^X+y^is calculated? it is then a random vector from 
the population N(y^,E^). This process is repeated for each pixel 
of the base classification and the random vectors thus generated 
are stored appropriately, i.e., so as to correspond to their 
simulated spatial location. 

The program requires as an input a classification map stored 
on a results tape. The results tape has the class statistics 


for p-dimensions also stored on it. The program, then, uses 
the results map and the stored statistics to generate a p-dimen- 
sional data set, which is stored on a user specified output tape 
in LARSYS format. 

Appendix I provides a mathematical derivation related to the 
generation of normally distributed samples. Appendix II provides 
the Fortran program listing for the simulation program. Appendix 
III provides the C program listing for the same program. 
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APPENDIX I 


Let and U 2 be two random variables independent and iden 
tically distributed Uniform (0,1). 

L 

Then let • (-2 In U^P cos 2 ttU 2 

k 

and “ (“2 In U^)^ sin 2 ttU 2 

then Z^ and Z 2 are independent and identically distributed 
normal (0,1) . 

Proof: 


f (U x ,u 2 ) 


1 0<Ui<l, 0<U 2 <1 

0 otherwise 


is the probability density function of two independent uniforms. 


* exp [-^(Z^ + Z 2 2 )] 



The Jacobian of the transformation is: 
J = - expOMZ.^ + z 2 2 )] 
f(Z x ,Z 2 ) = f(u x ,u 2 ) • I J I 

= exp [- *< ( Z 1 2 + Z 2 2 )] 

= 0 otherwise 

f(Z x ) * N ( 0 , 1 ) f(Z 2 ) ^ N (0 , 1 ) 


0 < [exp — *5 { Z ^ 2 + Z 2 2 ) ] < 1 


0 < ■= — arctan < 1 

2tt 



The side conditions give - ® < Z^ < ®, - ® < Zj < ®. 

Strictly speaking, Z^ cannot equal zero; however, prob (Z^ = 0) = 0 
as we are working with continuous densities. 

To test the effectiveness of the pseudo random vectors in 
the multivariate case, random vectors distributed N(0,I p ) were 
generated and then tested with a Kolmogorov-Smirnov test. Since 
the multivariate normal cdf is difficult to evaluate, the sum 
of squares was calculated and compared to the Xp 2 distribution. 

For sample sizes greater than 100, the pseudo random vectors 
were distributed properly. For sample sizes less than 100, the 
K-S test is not valid. Since we would generally (over an entire 
area) be working with more than 100 points per class, this was 
not pursued further. 

In addition, the sample covariance matrices wore tested for 
homogeneity against the true class statistics. For sample runs 
of up to 2000 points, there were not significant differences 
at the a = 0.10 level. 
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I DR EC I 7 1 * FLGT 

idrec! nii6i 3 * date;iii 

CONTINUE 

IDREC : 20 1 ■ NOLINE 
DO 145 II* 1. NOCHAN 
INEW * FETVC3I I I • 

WoftMH.IlW FRO 


SMRC226C 
SWRC2250 
SWRC226C 
SWRC227C 
SWRC228C 
SWRC229C 
SWRC230C 
SWRC231C 
SWRC232C 
SWRC233C 
SWK 02 340 


145 CONTINUE 

LIP * NOCHAN 


frocal: 1 1 2 . 1 NE W i 


00 150 II - LIP.NOLD 
DO 150 112 ■ 1.5 
FR3CAL:iI2.m « ' 


SWRC235C 

SWWC236C 

SWRG237C 
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me* SMITE FORTRAN 


PURDUE / LARS SCSI 


WRITTEN DTt DILI PFAFF .. , 

eoited rt* narman nuasher juke i^.irbc 







THIS PROGRAM GENERATES SIMULATED DAT-A RASED ON A 

SWiIHii‘?I8Sf5l! e i%«oS'r!8o:: , ' c 'r.SS'6. s t ..*?5 u t!5Ii .-f 
Tf'rSoSS^itsJliEBV^SMs cHostn *s . 


EAEC FILE ON TOUR DISK: 


GETDISK LARSTS 

UlOML TXTlIb CMSL 1 B FORTRAN SSP 37 C 

Eitfgtf \rAV*\w 

fili$Llle ? £actf e «Wf Sllti ^a^Sl 1 !^ peso 

RANOU nrtmtk 
START sprite 

THE PROGRAM RILL ASH FQR INFOAHATIC*. <UCH AS 
TAPE NUMBERS. MU NUMBERS, ..ETC. F«-“ Hfe.e UK, IT 
SHOULD BE EAST TO FCLLOr. „ 


SRSOGjIC 

sarcoczo 

SRRC033C 


SRR0009C 

SRRCOIOC 

SrRuOIIC 

SMRCQUC 

SmRCOUG 


variables USED IN TPRINT 


AREANO* 
B « 
DATA . 
DATVAL* 
ICAL « 
I OR EC * 
ISTART. 
LOGDAT. 
NOCHAN* 
NOCLAS* 
NOFLOS* 
NOPOOL* 
PNTCLS* 
l 


COVARIANCE STORAGE FCR FAC TCP I KG 

AREA NUMBER OF CLASSIFICATION 

COVARIANCE STCkAGE FOR MULTIPLICATION 

DATA POINT STORAGE 

LINE NUMBER ANC ROLL PARAMETEP 

CALIBRATION INFORM * [J ION „„ . , c 

IDENTIFICATION PECJRD STORAGE 

STARTING POINTS FO« GAUSS 

DATA POINTS IN LOGICAL FORMAT 

NUMBER OF CHANNELS IN CLASSIFICATION 

NUMBER OF CLASSES IN ORIGINAL STATISTICS 

NUMBER OF TEST MELDS 

NUMBER OF POOLEO CLASSES 

CLASSIFICATIONS ARRAV 

STATISTICS STORAGE 


INITIALIZATION 


IouIvIlcncc** V?$? i 1 »?: imd*t,logdati.:ical,icaii,:in-rt,ilini 
EQU IV ALENCE l FRO^ ALl 1 . 1 { , IDREC I 51*1 
OATA EOS, S. AM /'EOS »»1.C»C.C / 

DATA TES.Nfi, THREE /’TES • *NC '.'}•/ 


APPENDIX III 
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/ft******************************************** 

* 

* Read data from LARS Results Tare 

* And simulate data from it 

* usins the Box Muller relationship 

* 

#**«**•«*«*«««****«**««***#*****•*«***»#***«*** 

* 

* Swrite c has been translated from 

* the Lars Fortran Version of Suirite into 

* the lansuase 'o' for the Unix 0. S. 

* run on the DEC PDP-11/43 

* 

*************************************************** 

* Variables used in Swrite: 

* 

* a --- Covariance storage for factor ins 

* b == Covariance storage for multiplication 

* data -- Data point storage 

* nochan -= Number of channels in Classification 

* noclas Number of classes in original statistics 

* nopool -- Number of pooled classes 

* pntcls Classifications array 

* 2 ^ Statistics storage 

* 

* 

* compile with cc swrite. c -Ip -Ip /usr/1 ib/Pdsl ib 

* 

* 

* Initialise all variables used in swrite 

* External variables are available to all functions 

* within which thev arc declared 

*/ 

mai n< ) < 

extern irit noclas. no chan, norool, fdl, nonnts, nol i ne» recrium; 

extern int ier, fd2, 

extern float ens, uprcrCSl. lowerCSl; 

int ,i. ki fetvc3C53, debus. irifoC17J. rntcl si 10003, -i >, ipol; 

int i i , i done, iv» ma, mb, me; 

int ncch, nocomp, ictor, iond, ix, ibeg, 

i nt 1 nwrt, 12, i count, i r, n, intdat, i s tat 16 3. ip, io, it, nine, no, in; 
int nosam, ros, pf dl; 
int fd3, fd4, err5, 

char bufC 15003, *del im, *cl, *c2, pc 3, datoutt21003, obuft21003; 

char rnameC303, ttlMl 3, 

char mvfileC303, cmfileC303, opnamet303; 

float 2C6103, repnt, r; 

float bt63163; 

double rmeant31 3 16-3; 

float da ta 1 1 3 3 , 2214.103; 

float clapntC353; 

float mobufC303> 

double ivarC31 3 C63 C63, revar, imeanCSl 3 163, tl f 1 oat, t2f 1 oat, remean, semean; 
doubl e tempi. temp2. tempo, aC203. a2C63; 
double sart( ) , 1 og< ) , cos < ) , fmodi ) , s, t, f ; 

/« begin main program 

* if debug is set to minus one the following is 
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* printed out for verification on usr terminal: 

• nochan. nopool. f#tvc3. info, rntcls. nopnts. nolint. 2 . upur. lower. 
#/ 

debus <• li 

art ■ . OOOOl < 

/* read records 2. 4. 5. Set. off results tare*/ 

/* skip racords 1.3 */ 

readl to5( buf. f etvc3. 2 . info); 
noch = < < < nochan+1 )/2)*2); 
nocomp = ( noehan*< nochan+1 )/2); 
istop ~ nocomp * nopool; 
iend = nochan * nopool + istop; 
nosam = 4 * < < nopnts + 9)/4); 
if(nochan>5) £ 

pri ntf< "Number of channels is Xd but i ntarnal nochan) ; 

printf<" storase only allows 5\n"); 

pri ntf < "Executi on terminated abnornally \n"); 

exi t< ) ; 

> 

i f ( nopool>29) £ 

pri ntf< "Number of pooled classes is Xd but i eternal nopool ) ; 
printf<" storase only allows 29\n">» 
pr i ntf < "Executi on terminated abnormally \n" ); 
ex i t ( ) ; 

> 

i f < nopnts>~1000) £ 

pri ntf ( "Number of points- per line is Xd but i nternal nopnts) . 

pr i n t f < " s t ora 3 e o n 1 y all ow s 1 OOO \ n " ) . 

pr i n t f ( " Ex e c u t i o n t erm i na ted a b norma 1 1 v V n " ) ; 

exi t< ) ; 

printf < "nochan~Xd. nopool~Xd\n". nochan. nopool); 
i f ( debus -1 ) £ 

for< k=C; k<nochan; ++k) 

pri ntf ( " f etvcSCXdl ~ Xd \n" . k. f etvc3C kl ) ; 
for(k=0, k<i end; k = k+10) < 
for(j = k; ,K~k+9; ++■)> 
pri ntf ( " X f " . sC i ] ) . 
pr i n t. f ( " \ n ” i , 

> 

printf < "XnFi eld si?c: ")• 

printf ("\n line Xd to Xd with interval Xd". 

infoC41, infoCSl. infot&l); 

printf < "\n cols Xd to Xd with interval XdVn". 

infoC7], infoCS], .nfoC93). 

printf ( "Number of lines classified is Xdvn". r.ol i ne) ; 

pri ntf < "Number of points classified per line is Xd\n". nopnts) ; 

for< j=0; j<nochan; ++■)) £ 

printf < "upperCXdl - Xf " . i, urperC.il); 
printf <" lowerCXdl ~ Xf \n"» lowert.il); 

1 

read6rec< buf . pntols); 
i f < debus == -1 ) £ 

printf < “XnFirst line of r ecord Xd f ol 1 ows\n“. reenum) ; 
f or < J=0; i<=nopnts; j = i + 10) < 
for< k=J; k<= i+9; ++k) 

printf<" Xd". pntclst kl ) ; 
pri ntf < "\n" ) ; 

> 


■i.J t 
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> 

/* m-i I; e a o. v of statistics Array 2 */ 
for< ix*l# ix<®i end; ++ix) 
z2Cix3 - 2fix-13; 

/* Create outrut file */' 

sets( "VnSreei fv lincmnUr outrut file pathname ", orname), 
fd2 ~ or eat( orname, 0777) . 
if ( fd2<0) 

pri ntf< “Cannot or cat lr outrut file ?.s\n", orname) ; 
pfdl - 1; 

scts( "Srccify simulated data PDS file ra thnamc: ",rname); 

3ets( "enter a dO char ti tic: \n“, ttl ) ; 

sets< " ' nSreci fy mean vector file pathname: ",mvfile)> 

sets( "VnSpooi fy covariance matrix file pathname: ”,cmfile); 

rdsoren( rf dl, rname); 

putfmtlrfdl, 1, nosam, nolino.S, nochan); 

rutttl ( pfdl, ttl); 

fd3 - or ea t < mv file, 0777 ) ; 

i f < f d3<0 ) 

rrintf < “Cannot creat mvfilo Xs (MAIN) Vn", mvf i 1 e) ; 
f dd ~ cr ea t < cm f 1 1 c • ' *'777 ) ; 

1 f < f dd<C») 

pr i ntf( "Cannot creat cmtile %s (MAIN) vn", cmf i 1 e) ; 

/* write output pa sc one #/ 

c 1 ~ " ♦Da ta Simulati on U s ins B ox-M ullcr Relationship* 

rrintf(fd2, "\p\n\n\h%s\n%s\n%s\n“, dclim, cl, dclim), 
pri ntf ( f d2, "\n Line Xd to line Xd with interval Xd \n“, 
infold], i nf ot*53, i nf oC*l ) ; 

pri ntf ( f d2, "\n Column X<i to column X d with interval Xd vn", 
infoC7], infold], infot9]); 
rri ntf < f d2, "\ri channels used vn"), 

f or< ix=l» ix<~noch*n; **ix) 

pri ntf ( f d2, " Xd Xf - Xf\n", f etvc3tix-ll, lowert ix-13, 
urr er C i x- 13); 

/* output new rase character */ 
pri ntf < fd2, “\01d"); 

/* Factor covariance matrices */ 
ibes =■ 1; 

f or ( ix~l ; ixOnor ool ; ++ix) C 
idone ~ ibes ♦ nocomr- - 1; 
k - 0; 

f or ( i y=i bed; iy<~i done; **i y) t 
k ~ N + 1; 
alk] = sCiv-l]; 

> 

/♦ call function to perform factor i ns */ 
mfsd(a); 
i f ( i er == - 1 ) 

pri ntf < fd2, “Error -l\n“); 
i f ( i er>= 1 ) 

pri ntf <fd2, "Error st l\n"!; 
i f ( ( i er == -1 ) II tier >= 1>) 
exi t< 1 ); 
k = 0, 

f or( i v=i bes; iv<=idone; ++iv) t 
k » k + 1) 

2 liy-l] = aCkl; 

> 

ibe9 = ibe» ♦ nocomp; 
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> 

/* initialize random number senerator */ 
srand* 1 ) ; 

/* initialize arrays */ 
for*ma E l; ma<-noclas; ++ma) < 
clepntCmel = 0. 0; 
for*mb-l; mb<=nochan; ++mb) C 
imeanCma 1 Cmbl = 0 0; 
rmeantmal Cmbl - 0 . 0 ; 
f or<mc=l; mc<~nochen; ++mc) < 
ivarCmelCmbKmcl - 0 0; 

> 

> 

> 

t = 25. 0; 

lnwrt - 0; 

/* this while loot 1 ic repeated tor each line in the classification */ 
while*recnum =•* 6) < 

1 nwr t = 1 nwr t + 1 • 
s = lnwrt, 
f * fmod* s, t ) ; 
i f < f =-- 0 O ) 

rri ntt'( ">.d Lines out of Xd are c omrl eledvn", 

1 mwrt, nol i ne) » 

12 - 0; 
i count ~ 4 

f or ( i x~ 1 ; i x<- n o r n t s ; ++ 1 x: > C 
i count - i count + 1; 
i2 - rntcl sf i: -1 3; 
irol - <i2-l) * nochun, 
l tea “ (i2-l> * no com r-; 

I-: = ibes; 

for* iv*l; iv<~noch4n; ++iv> < 
for* i 2~1; 1 2<1* i v; t+is) < 
k - k + 1, 
bCivlCirl - sCL-ll; 
if* i y •- i 3 ) 

Mi:H i y3 - 0 0; 


for*iY~l» iv<- noch; ++i r) < 
e2Civl ~ r^ndU, 
a C i v 3 - rend*); 

■i2 C i y 3 * a2C i v 3 / 32767 , 
ali/3 ~ ativl / 32767 • 

ativl = strt* -2 0 « 1 os* a2C i y3 ) ) *■ eos(6. 28318 * aCix]); 

> 

clapntt 123 - clapntt i21 + 1.0; 
for * i y- 1 » i Y<~noche n; ++iv) { 
da tat i y 3 ~ 0. 0, 

i = nopool * nocomp + ir ol + iv; 
for* i z=l , i 2 <~nochari, ++i 2 ) 

datatirl - datativl + bCivltizl * atizl; 
dataCiY] = dataCirl + 2tie— 13; 
intdat = dataCiv] + 5, 

l f ( i ntdat<0) intdat = 0, 
i f < i ntdat; 255) intdat ^ 255; 
istattivl = intdat, 
pos = <iY-l)*nosam + icount - 5> 

if(POS>2100) printf*"datout internal buffer full\n M ); 


19 




dotoutt po*3 ■ lntdot tc 0377; 
fort lz-1; iz<-6; 2 ) 

d4tout[pos+lz3 - 0; 

> 

f or< 1 1-1; 1 l<-nooh4n; ■♦■■*■1 1 ) < 

imtonC 123 C 1 1 3 ■ ImoonC 123C til ♦ i*t4tC113; 
fort JJ-ii; J JOnoohen; ♦♦j j) < 
ttflirO ■ l*t4tC113; 
t«*IPl * i*t4ltjJ3) 

lvorC 123 C 1 1 3 C Jj 3 ■ 1 v*r C 123 C 1 1 ] C j J ] + t«mpO * tempi; 

> 

> 

> 

il-O; 

fort 1 y« 0; lY<no*4m; ++1 y) < 
f ortiz-0 ; i2<nooh4n; ♦■*■1 j) 

obufCii++3 - dotoutC nos4m*i 2 ♦ 1y3; 

> 

Pd*lPo*< pfdl. lnwrt-1); 

i f t putl 1 n*( pf dl> 1 nwr t-l» obuf. no*4m*noch4n ) !- O) 
printft "Putline to PDS output f4ll#d\n"); 
r«4d6r«ct buf. pntcl*); 

> 

/* *nd of while loop */ 

/* 

fort 12=1; i2<=18; -m-12) < 

printft fd2, "\n %d \n", 12 ); 
fort i 1 = 1; 1 i<=-noch*n; 1 ) C 
f ort j j=il; j ]<=noch4n; j j ) < 

pri ntf t f d2. ” % f ", iv4rCi23Cii3t jJ3>; 

> 

pri ntf ( f d2< "\n" ); 

> 

> 

fort 12=1 ; i2<“18; i 2 > < 
pri ntf ( fd2, "\n"); 
f or ( i i-1; 1 K-nochan; ++i i ) < 

pri ntf t f d2< " */.f \n", imoonC i23 C i i 3 > ; 

> 

> 

*/ 

fort ip^l; ip<=noclis. ♦♦Ip) < 
fort io=l; lo<=nochin; ++io) < 
i ft clopntC ir]>0. 0) < 

1 1 f 1 04 1 = ime&nCip3Clo3; 

t2fl04t = cl4rntCip]; 

rm«4nt 1 p 3C i o] = tl f 1 o* 1/ V2f 1 oa t; 

> 

f or t i t=i o; i t<=noohan; ++i t) < 

1 f t clipntt 1 p3>1 0) < 
repnt * clarntClpl; 
r«v*r = ivarC 1 p 3 C i o3£ i t3; 
rtnun - imeanC 1 p 3C 1 o3; 
semean - imeanC i p3C i t3» 
t*mrO = re<J4r/tr«pnt-l. 0); 
tempi = r«msan/r«pnt; 
t«mp2 = s«m«an/trepnt-l. 0); 

Wart 1 p] C 1 o3 C i t3 - t«mpO - ( tempi* t«mr2) . 

WarC 1 p 3C 1 13£ 1 o3 - WarC 1 p 3C i o] C 1 t3; 

> 
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> 

> 

> 

/* output results */ 

for(ip»li ip<«nool4*i ♦♦I p) < 

printf(fd2. " Class number Xd Xf poi nts\n\n\n". 
i p» clarntt i r3 ) ; 

printf < fd2, " Aotual Simul*ted\n" ) 

printf<fd2. " Mean m«*n\n"), 

for< ix~l; ix<=noch*ni -*"*-ix) < 

nine * nocomp * noclas ♦ <ip - 1) * nochan. 

pri ntf ( f d2, " Channel Xd < Y.6. 3f - X6. 3f ) X6 3f Xfe 3f\n", 
f etvo3C ix-1 3. lowortix-t3. urrertix-13, z2Cninc+ix3. 
rmeantiP3Ci>;3)» 
mobuftix-13 ~ rmcanC i r3C ix3; 

> 

errS = wri te( fd3, mobuf , 4*nochan>; 
if < errSCO) 

printf ( "Error- occurred writinz Mean Vector file (MAIN) \n")» 
pri ntf < f d2. "\n\n\n\n\n Actual Covariance MatrixVn"); 
fori no=l; noOnocomr. ++no) < 
nine = <ip-l) * nocomr; 
aCno3 - z2C m nc+no3. 

3 

wr tm tx < a, f c tv c3 ) < 

pri ntf ( f d2, n' n\n\n Simulated Covariance MatnxVn"); 
no * 0. 

foriio^li i o<“nochan. ++i o) C 
f or ( i n^-1 i i n<~i o. ++i n) < 
no — no + 1; 

at no 3 - ivartir 3Cio3t in3; 
mobuf t no-1 3 ~ i"arC j r 3 1 i o3 1 1 n3; 


> 

wrtmtx(a. fetvc3)> 

err5 ~ wri tc( fd4. mobuf, 4ttriocomr) . 

i f < err 5 <0 ) 

printf < "Error- occurred writms cm file (MAIN) \n“); 
printf < f d2< "\014" ) • 

> 

/* cexit terminates activities on or on files and flushes 
* the output butter cexit is part of the c library */ 
printf ( "swri to. c is t i ni shedvn" ) ; 

printf ( "Vn The simulated data (in PDS format) is at Xs vn",pname>. 
printf("The linernntor output file is at Xs vn".orname); 

Printf ("The Mean vector file is at Xs \n“< mvf 1 1 e) . 
printf("The Covariance Matrix file is at Xs \n".cmfile). 
pdsolose(pfdl)> 
o«xi t( ) . 

> 

/* end of main program */ 

/**»##»***#t«**it**e*****e*«*f*f*«tt*f***»**»*»»»*»«**«*/ 

/* function to read records 1 thru 5 follows */ 

r eadl to3( buf . f etvc3, word, info) 
int *f etvc3# *i nf 0 / 
char *buf; 
float *word; 

C 

extern int noolas. nochan. norool, fdl. nornts. noline; 
extern float uprerC33. 1 ower[53; 


■ Vi 


. '* l L pAo,, 


IS 
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Int errl, J, k, reenum, stetsi *e; 

oh *r lbufC41> 

floet ibmdect ) , 

fdl - oren( "/deWrmtO", 0>; 

i f ( f dl<0) 

printf ( "Cannot open 9 trk r.ere flit (READ1T05) \n">; 

/* skip reoord 1 *f 

errl * re*d< fdl. buf, 1300), 

1 f ( trr 1 »■ -1) 

rri ntf< "Error occurred in reading record 1 (READ1T03) \n")i 
/* reed record 2 */ 

errl * reed< fdl, buf, 1300); 
if (errl *« -1) 

printf ("Error oocurred readins record 2 (READ1T05) \n"); 
nodes * bufC193; 
nochan - bufC233, 
nopool = bufC313, 

f or < J*>0, k*-33; .Knochan; k-k+4) 

fetoo3C33 “ buf C 1:3; 

f or< k*32+4*noohan> j~0; j<2Pnochan; ♦♦.i, k-l«+4) { 
tbuf CO] * buf CIO) 
tbufCl] - buf C k+1 3 ; 
tbuf C23 = buf C l;+23; 
tbuf C 33 - bufCM-33; 
i f < Knochan) 

lower!)] ~ i bmdec( tbuf > . 
else uprerC f-nochan] - i bmdec'. tbut ) ; 

> 

/* skip record 3 */ 

errl =■ read( fdl , buf , 1500) , 
l f < err 1 — - 1 ) 

pr i ntf< "Error- occurr ed readina record 3 (READ1T05) \n“); 
reenum - bufCllJ. 
whi 1 e< r eonum -- 3) i 

errl - read', fdl > buf , 1500) , 
l f < errl --- -1) 

rri ntf( "Error occurred reading record three < READ1 T03 ) \n" ) 
reenum - buff 11]. 


/* Handle record 4 

statsise “ nochanc norool + ( noeharipf nochan+1 ) /2 ) enorool i 
for<k“0, )“16, IXstatsi se. ++k, j~i+4) < 
tbuf CO] ~ buf C O, 


> 


tbuf C 1 ] " bufCi-MJ. 
tbuf C2] - bufC)+23. 
tbuf C3] » buf C )♦:?•], 
wordCk] ” l bmdec'. tbuf ) , 
> 

errl * read( fdl , buf 1300) ; 


if (errl -1) 

printf ( "Error occurred readina record 5 (READ1T03) \n"), 


nopnts * 
noline *■ 
infoC4] * 
infoC53 “■ 
infoC63 - 
infoC73 * 
infoCS] * 
infoC9] 


( buf C 183<':8 I 
< buf C 22 3 <<8 I 
( buf C383 < :8 I 
( but C423< 18 I 
( bur C463«8 I 
( buf C30J<<8 I 
< buf C34 3< ;8 I 
(bufC383C'8 I 


< buf C 193 

< buf C233 

( buf C393 

< buf C433 

< buf C473 
(bufcsn 
( b u f C 35 3 
( buf C593 


t, 0377 ) ) , 
& 0377)); 
t, 0377)) 
0377) ) 
0377) ) 
0377) ) 
0377) ) 
0377) > 


lnt nool**, norool. noohan, fdl. nornt*. noline; 
float urr*r£53, lower £53; 
float ib*d*c(buf> 
ohar *buf; 

< 

int k. l,m, tint; 

char nbuf £23. to her temr3> bvt*£43, sl«n; 
f 1 04 1 word; 
m * 0; 

*i*n - 0000; 
i f < buf £03 < 0) 
fi*n ^0200; 

ti nt*< < < buf £03 *< 0177) - 64) * 4); 
whl 1 e< buf £ 1 3 :> 0) < 

buf £ 1 3“< < buf £ 1 J«1 ) b 0376); 

♦+m; 

) 

It * 8-m; 
nbuf £03“buf £23; 
nbuf£ 1 3 •'buf £33; 
f or< 1 "t 1 ; 1<-It; ) 

nbuf C03"< ( nbuf t ! 3»1) b 0177); 
nbuf E 1 3~< < nbuf £ l 3»1 ) b 0177); 

> 

f or< 1-1 , l<~m; ) < 

huf£23-< < buf E23CI ) b 0376); 
buf £33"< < buf £ 33<<1 > b 0376); 

> 

buf T 1 3~< buff 1 3 I nbuf £03); 
buf £23-'. buf £23 I nbuftlS); 
ti-)t - ti nt-m+128; 
ifitint < 0) tchar - OoOO, 

*l*e iff tint > 235) tchar - 0377; 

site tchar - tint; 

buf £ 1 3'( buf £ 1 3 b 0177), 

t em r 3 - tchar«7; 

tefiir-3 ~ tomr3 b 0200; 

buf£13 - temr3 I buftl3> 

buf£03 “ tehaK'M, 

buf £03 - buf £03 b 0177; 

buf £03 - buf £03 I nsn; 

bv tc£03 “ buf £ 1 3 b 0377, 

bv te£ 1 3 buf £03 b 0377, 

bvt*£23 - buf £33 b 0377; 

bvteE 33 ■* buf £23 !< 0377, 

racM bvte, bword); 

returnfuiord) ; 

> 

/* the function rad- taker the 4 8-bit character bvter 

* that h*M* been rearransed bv ibmdec and racl-r 

* them into a floating roint word */ 
r*ck(brte. cword) 

char *byt«, *cword; < 
int i; 

for< j"0; i<4; ■*■+)> 

cword£i3 ’* bvte£ ;3, 

> 

/■*##**##******#*********-**********r*******************c 

/* function t<* r <"i ,1 r ••cord r. tel loot *■/ 
read6r«c< buf, rntcl:) 
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lnt eentolti 
ch*r ebufi < 

extern lnt norntt, recnirv i dlt 
int i> n. errl, k > 
n • 0* 

errl»reed< #dl. buf, 1500) ; 
i f < err 1 »» -1) 

wintfl "Error occurred reedlne record b\n“U 
rtcmmabufUDt 

#or< J"20, k«Oj k<nornlsj k*+. C 

rntcliCnD^bufCJe-lJ & 0377 ) j 
♦♦iw 
> 

> 

lnt r»cnwi 

f* Mftd feotort th« Mlrloit recced bv th« iMin rro*T4*. 

• Mftd It * direct trencletlon from the Lers 

• fortren version of the seme neiee */ 

»f td<*> 

double *e* < 
extern lnt noohen, ler. 
extern float erf; 

lnt kriv, k, lnd, lend, 1, lenf, 1. llnd; 

double tol, dtuei, drlv, y, feb£< ); 

double done, ter t< ), x; 

done * 1 0; 

ler * O, 

kriv * 0; 

i or< k"i; kC-nochan; ♦♦): ) < 
kriv r kriveli, 
lnd ■» kriv. 
lend » k-1; 
y * ert«e( kelv), 

lol «■ febc<Y). 

for< i"k; K-nochan; ♦♦! > < 
dime *■ 0 0, 
if< lend 1 - 0) < 

f or< 1-1, K-i e nd, eel ) <. 

1 4 n f - k r i v- 1 , 
llnd ind-1; 

d £ urn - ( dium*<at lenf 3«el 1 1 nd3 ) ) ; 

> 

> 

diu<n - at 1 nd3-dcu»; 
if< < i-k) •« 0» < 

i f < < ( dsue-tol ) <" 0.0) kk <dtum > 0.0) kk (ier<»0)) 
ler • k-1; 

1 f< ( < dfuei-t . O <- 0 0) kk (dtum <■ 0 0) ) C 
ler * -1; 

rri ntf < "E r ,•> -lvn"»; 
exlt(l); 

> 

x • dsuei; 
driv m ter t<x>; 
etknvl ■ dr»v> 
driv ■ done/del v; 

> 

1 f < < i-k) •- 0) 

eCindl " dtu*«driv; 
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ihd ■ ind+ii 

> 

> 

> 

int 4 <tr> 

f i04t •fS) 

/*«M*»*******i>«***f*#**m*»#tm*«#**##*#**## M ***t«/ 

/* the write *iUir functl n <riU: to the outrun 

* f41< th< 1 owin' half vf t. ntchinf nochifi covorience mVtix 

* which is r*rscd through the NrM«i«r 4 . #/ 

wrU»tK<4. fetveS) 

lot *fetvc3, 
double 04 . c 
extern int fd2, nochin, 
oxtern float urrerCSl, 1 owertSJ. 

4 nt i, k» m> 
m "■ 0; 

rrintf (fd2»"'n 'sreotril ")> 

fori i’-O > -Krtoehan; *♦ . 1 ) 

rrintf ( f 42. "\7 . 3i - ”, 1 odct t i 3 ) , 

rr intf < f <12, ” vn Bind ”), 

t or ij“ O; > 'nochifi. ♦♦ 1 ) 

rrintf<fd2, "\7 r.f ”, urrc* l 4l>, 

for Knochan. ♦♦.<) < 

r» »nlf< fd2, n' r. \? ?>( -vn". 1 owert J 1 > . 
rrintf<fd2>" %7.*>f ",urrerljJ>, 

for<l— 0; l;<- », ♦♦!,) <. 

rrintf < fd2- "\7 3f ".aH"*-*!); 

> 

> 

) 

int fdi; 

/ , ***********<‘*Mt*o*(>*»»»^«*#»»t*«#*(HHl**»*(nn l **( M nK ( ( ((t / 

»etf < r, r ) 
ch4 r *r, *n < 
ehor c, *rli 
rrintf(r), 
el * s; 

while* <c » seteharO) ‘w -\n') 

*rl++ ” ci 

*rl » 0; 


